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We investigate higher than the first order gravitational perturbations in the Newman- Penrose 
formalism. Equations for the Weyl scalar ipi, representing outgoing gravitational radiation, can be 
uncoupled into a single wave equation to any perturbative order. For second order perturbations 
, about a Kerr black hole, we prove the existence of a first and second order gauge (coordinates) and 

■ tetrad invariant waveform, tj)! , by explicit construction. This waveform is formed by the second order 

piece of i/'4 plus a term, quadratic in first order perturbations, chosen to make 4'i totally invariant 
and to have the appropriate behavior in an asymptotically fiat gauge, ijji fulfills a single wave 
' equation of the form Tipi = S, where T is the same wave operator as for first order perturbations 

and 5* is a source term build up out of (known to this level) first order perturbations. We discuss 
the issues of imposition of initial data to this equation, computation of the energy and momentum 
radiated and wave extraction for direct comparison with full numerical approaches to solve Einstein 
equations. 
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I. MOTIVATIONS AND OVERVIEW 



The prediction of accurate waveforms generated during the final orbital stage of binary black holes has become 
a worldwide research topic in general relativity during this decade. The main reason is that these catastrophic 
astrophysical events, considered one of the strongest sources of gravitational radiation in the universe, are potentially 
observable by LIGO, VIRGO and other interferometric detectors. For its strong nonlinear features this black hole 
merger problem is only fully tractable by direct numerical integration (with supercomputers) of Einstein equations. 
^ ' Several difficulties remain to be solved in this approach such as the presence of early instabilities in the codes 
for numerical evolution of Einstein theory and finding a new prescription for astrophysically realistic initial 
^ data representing orbiting black holes 1^,^. Meanwhile, perturbation theory has shown not only to be the main 
approximation scheme for computation of gravitational radiation, but also a useful tool to provide benchmarks for full 
numerical simulations. From the theoretical point of view perhaps the more relevant contribution during the nineties 
in perturbative theory has been the "close limit approximation" Q|. It considers the final merger state of two black 
holes as described by a single perturbed one. This idea was applied to the head-on collision of two black holes and the 
' emitted gravitational radiation was computed by means of the techniques used in first order perturbation theory around 
a Schwarzschild black hole. When the results of this computation have been compared with those of the full numerical 
integration of Einstein equations the agreement was so good that it was disturbing g. This encouraged the significant 
effort invested into the development of a second order Zerilli formalism of metric perturbations about the Schwarzschild 
background. The method was successfully implemented with particular emphasis on the comparison with the fully 
numerically generated results. In the case of two initially stationary black holes (Misner data) the agreement of the 
results is striking |6| . Second order perturbation theory confirmed the success of the close limit approximation with 
an impressive agreement in both waveforms and energy radiated against the full numerical simulations. There has 
been a tantamount success in the extension of these studies to the case of initially moving towards each other black 
holes 0, and for slowly rotating ones ||] (See Ref. for a comprehensive review). 

All the above close limit computations are based on the Zerilli jlo| approach to metric perturbations of a 
Schwarzschild, i.e. nonrotating, black hole. This method uses the Regge- Wheeler [|ll| decomposition of the met- 
ric perturbations into multipoles (tensor harmonics). Einstein equations in the Regge- Wheeler gauge reduce to two 
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single wave equations for the even and odd parity modes of the gravitational perturbations. There is, however, the 
strong belief that binary black holes in a realistic astrophysical scenario merge together into a single, highly rotating, 
black hole. There is also concrete observational evidence of accreting black holes that places the rotation pa- 
rameter as high as a/M ~ 0.95. Finally, highly rotating black holes provide a new scenario to compare perturbative 
theory with full numerical integrations of Einstein equations. 

The Regge- Wheeler- Zerilli techniques cannot be extended to study perturbations on a Kerr black hole background 
(see Ref. B for the slowly rotating case). In this case there is not a multipole decomposition of metric perturbations 
(in the time domain) and Einstein equations cannot be uncoupled into wave equations. A reformulation of the 
gravitational field equations due to Newman and Penrose , based on the Einstein equations and Bianchi identities 
projected along a null tetrad, allowed Teukolsky fl^ to write down a single master wave equation for the perturbations 
of the Kerr metric in terms of the Weyl scalars il)^ or i/>o. This formulation has several advantages: i) It is a first 
order gauge invariant description, ii) It does not rely on any frequency or multipole decomposition, iii) The Weyl 
scalars are objects defined in the full nonlinear theory and a one parameter perturbative expansion of it was proved 
to provide a reliable account of the problem . In addition, the Newman-Penrose formulation constitutes a simpler 
and more elegant framework to organize higher order perturbation schemes as we will see in the next section. 

Since the seventies the Teukolsky equation for the first order perturbations around a rotating black hole has been 
Fourier transformed and integrated in the frequency domain for a variety of situations where initial data played no 
role (see Ref. for a review). Very recently it was proved [|l^,0 that nothing is intrinsically wrong with the 
Teukolsky equation when sources extend to infinity and that a regularization method produces sensible results. In 
order to incorporate initial data and have a notable computational efficiency, concrete progress has been made recently 
to complete a computational framework that allows to integrate the Teukolsky equation in the time domain: First, 
an evolution code for integration of the Teukolsky wave equation is now available and successfully tested [ pO| . 
Second, non conformally fiat Cauchy data, compatible with Boyer-Lindquist slices of the Kerr geometry, began to be 



studied with a Kerr-Schild [£1)22| or an axially syinmetric |23| , |24| ] ansatz. Finally, an expression connecting ip^ to only 
Cauchy data has been worked out explicitly |2^j2^,|2^ . 

Assuming that we can solve for the first order perturbations problem, we decided to go one step forward in setting 
the formalism for the second order perturbations. As motivations for this work we can cite the spectacular results 
presented in Ref. |^ for the head-on collision and the hope to obtain similar agreement for the orbital binary black 
hole case in the close limit. Second order perturbations of the Kerr metric may even play a more important role in 
this case since we expect the perturbative parameter to be linear in the separation of the holes ]2^ ] while in the head 
on case it is quadratic in the separation . The nonrotating limit of our approach will also provide an independent 
test and clarify some aspects of Ref. Q results. High precision comparison with full numerical integration of Einstein 
equations using perturbative theory as benchmarks is also one of the main goals in this program as well as a the 
development of a tool to explore a complementary region of the parameter space to that reachable by full numerical 
methods. An important application of second order perturbations is to provide error bars. It is well known that 
linearized perturbation theory does not provide, in itself, any indication on how good the perturbative approximation 
is. In fact, it is in general very difficult to estimate the errors involved in replacing an exact solution of the full Einstein 
equations with an approximate (perturbative) solution, i.e., to determine how small a perturbative parameter e must be 
in order that the approximate solution have sufficient accuracy. Moreover, first order perturbation theory can be very 
sensitive to the choice of parametrization, i. e. different choices of the perturbative parameter can affect the accuracy 
of the linearized approximation [p8| . The only reliable procedure to resolve the error and/or parameter arbitrariness 
is to carry out computations of the radiated waveforms and energy to second order in the expansion parameter. The 
ratio of second order corrections to the linear order results constitutes the only direct and systematically independent 
measure of the goodness of the perturbation results. 

In the next section we extend to second (and higher) order the Teukolsky derivation of the equation that describes 
first order perturbations about a Kerr hole. To do so we consider the Newman-Penrose formulation of the Bianchi 
identities and Einstein equations, make a perturbative expansion of it, and decouple the equation that describes the 
evolution of second (and higher) order perturbations. This equation takes the following form 

f^P^^^ = S[,p^^\dt^l^^% (1) 

where 'i/;=(p(*'))^'*'04, is the same (zeroth order) wave operator that applies to first order perturbations (see Eq. 
(p^)) and S* is a source term quadratic in the first order perturbations (see Eqs. (P)-(pd])). 

In Section III. A we describe how to compute the source, appearing in Eq. (|l|), in terms of solutions of the wave 
equations for i)'"^^ or i/jq"'^'' only, which are the objects we directly obtain from the integration of the first order Teukolsky 

equation. Sec. III.B discusses the issue of building up ip)^ and dti^X out of initial data (that we assume are given to 
first and second order). In section III.C we recall the equations for the computation of the second order total radiated 
energy and momentum. 
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Higher than first order calculations are always characterized by an extraordinary complexity and a number of subtle, 
potentially confusing, gauge issues mainly due to the fact that a general second order gauge invariant formulation 
is not yet at hand in the literature. In general, gauge invariant quantities have an inherent physical meaning and 
they automatically lead to the simpler and direct interpretation of the results. In the Newman-Penrose formalism 
one has not only to look at gauge invariance (i. e. invariance under infinitesimal coordinates transformations), but 
also at invariance under tetrad rotations (see Sections IV. A and IV. B). More specifically, the problem here is that the 

waveform il)^^'' in Eq. (^) is neither first order coordinate gauge invariant nor tetrad invariant. The question that 
arises therefore is whether can be unambiguously compared with, for instance, full numerical computations of the 
covariant tp^""^. To handle this problem we build up a coordinate and tetrad invariant quantity up to second order, 
tpj- , which has the property of reducing to the linear part (in the second order perturbations of the metric) of '04 hi 
an asymptotically flat gauge at the "radiation zone" , far from the sources. This property ensures us direct comparison 
with i/il^"™ by constructing iIj^^'' + ip'"^]. In Sections IV. A. -IV. C. we give an explicit and general prescription for the 
construction of second order gauge and tetrad invariant objects representing outgoing radiation. To do so we impose 

(2) 

the waveform ipj to be invariant under a "combined" transformation of both the coordinates and the tetrad frame 
to first and second order. The resulting second order invariant waveform can then be built up out of the original 

(2) (2) 

■04 plus corrections (quadratic in the first order quantities) that cancel out the gauge and tetrad dependence of 04 . 
Finally, in Sec. V, along with a short summary, we discuss the astrophysical and numerical applications of our result. 
We end the paper with three appendices: Appendix A refers Sec. III. A and contains explicit formulas to compute the 
first order perturbative Newman-Penrose quantities (Weyl scalars, spin coefficients and perturbed tetrad) in terms of 
the first order metric perturbations needed to build up the source term in the wave equation for 0j . Appendix B 
refers to Sec. III.B and contains formulas to compute the second order spin coefficients in terms of the second order 

(2) 

metric perturbations and product of first order perturbations needed, for instance, to build up 0) in terms of initial 
data. Finally, in Appendix C we explicitly give the expressions to build up the gauge invariant waveform holding in 
the Schwarzschild limit case, i.e. for a = 0. 

Notation: In this paper we use Refs. [t30yi3[ conventions. Background quantities carry the (0) superindex if needed 
for clarity and are all explicitly given in the cited references, while superindices (1) and (2) mean pieces of exclusively 
first and second order respectively, for instance, we expand = 0;*^°^ -I- ■0^^-' -I- 0^ 
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II. DECOUPLED EQUATIONS FOR HIGHER ORDER GRAVITATIONAL PERTURBATIONS 

Let us consider the following two of the eight complex Bianchi identities written in the Newman-Penrose formalism 
(projected along a complex null tetrad) Ch. 1.8 (see also the Appendix A) 

{D + Ae- p) 04 - (5 + 47r + 2a) 03 + 3A0J2 (2) 
= 47r[(^ - 2r -I- 2a) T„^ - ( A -I- 27 - 27 -f 71) T^sts - A (T„; -I- Tmwi) + oT„„ + vT^] , 

((5 + 4/3 - t) 04 - (A + 4m + 27) 03 + 31/02 (3) 
= 47r[(^ - r + 2/3 + 2a) r„„ - (A + 27 + 2J1) T„^ + v (r„; H- T„-crn) + vTrnm - AT„m] 

and the following one out of the eighteen complex Ricci identities 

(A + M + 71 + 37 - 7) A - (J + 3a + ^ + TT - r) i/ + 04 = 0. (4) 

Here D = l^d^,, A = n^a^, 8 = m^'^^,. 

In what follows it is convenient to define the operators 

d3=(5-H3a-H;5-h47r-r) , ^4=^ (A + 4/i + p + 37 - 7) . (5) 

In order to find a decoupled equation for i/'4 we operate with|^ ^4"^ on Eq. (||), with dg"' on Eq. (||), and then 
subtract to obtain 



'^Here we use operators defined on the background instead of (^ for the sake of simplicity. 
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{D + Ae-p)- df (5 + 4/3 - t) 



■77(0) -7(0) ■ 



i>2 



= r [matter], 



'04 



4°^ (A + 4^ + 27) - (5 + 47r + 2a) 



V'3 



(6) 



where T [matter] is defined in Eq. ([lOl ) below. 

In the above equation ip^, ^3, v and A vanish on the background, i.e. on the Kerr geometry, but so far this equation 
is exact, no perturbative expansion has been made yet. Let us now think how to use Eq. (^) in a perturbative scheme. 
In this context, the superindex (p) appearing in the formulae below stands for a sum over all perturbative orders from 
p = 1 up to p = n — 1 (i. c. ) where n — 1, 2, ... is an arbitrary order we want to study. 

To fix ideas let us first discuss second order perturbations, n — 2. The procedure for higher order perturbations 
will be clearly analogous. We want to have an uncoupled equation for V'i^^- Since '04°^ = 0, the operator in the 
first bracket on the left hand side of Eq. (^) is needed to zeroth plus first order. The zeroth order acts on ■04^-' and 
generates the same wave operator as for the first order perturbations. The first order operator in the first bracket 
on the left hand side of Eq. ^ acts on ilif^ and its result can be considered as generating an additional source term 
since it is supposed we have already solved for the first order perturbation problem previously. The second bracket 
on the left hand side of Eq. (ph can be considered as a pure source term as well since its zeroth order vanishes 



dr(A + 4;. + 27f^ 



df (^ + 4^ + 2a)*°^ =0, 



(see Ref. '30| for an analogous proof) and then we have to consider "03^^ only to first perturbative order (in 
general, to all lower perturbative orders than the one considered). The last bracket on the left hand side of Eq. (H) 

includes terms depending on v^"^^ and A^^' since "^^^^ —M/{r — iacosi?)'^) , is non vanishing. To get rid of these 

second order spin coefficients we use Eq. (^ multiplied to the left by ■02"'' 



-df At") 



(d4-3M)'"~'^ A(rt 



where we have made use of the pure zeroth order relations At°)'02°) = — 3/it'')?02°'' and 5 °V2°'' = —iir'^^^ip^^ coming 
from the Bianchi identities. The above result allow us again to write the terms depending on -02 as source terms. 
We finally obtain the equation that describes the n-th order perturbations 



{df {D + U~ p)(°) - df (^ + 4/3 - r)(") - 3^°)} V4"' = 54[^("-''), + T[: 



T [matter! 



(7) 



where 



(8) 



and the source terms are (where brackets represent operators) 



54-E [4"^^ + 4/3-r)( 

- (A + 4a* + 27) 
+3[df i.("-P)-df At"-'')] 4^' 



-df {D + Ae-pf 



df (A + 4A* + 27)^"-''^-df 



47r + 2a\ 



p) 

[n-p) 



4'^ 



(9) 



(0) 



(d3 



o (p) 



(d4 - ip) 



{n-p) 



Xip) 



and 
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n-1 



0=1 



r[matter] = ^ (S - 2t + 2a) T^g - (A + 27 - 27 + p)^ 



T— 



+^3 



(A + 27 + 27I)("-^) T^rt _ (J _ T + 2^ + 2a 



("-p)y(p) 



(10) 



where T^^ = {Tf,ynt'Wi}')^P\T^^ = (T^^m^m^j^P) and Ti?'? = {T^un^'n'')^P\ Note that in our formaUsm we have 
taken into account matter terms in order to be used in future computations including an orbiting particle or an 
accretion disk around a Kerr hole. By summing up over all n— orders in Eq. ^ one should be able to recover 
solutions to the full Einstein Equations. 

Note also that if one wants to act on ip^'^'^=p~'^ip'^^^ rather than tp^^ one should rescale all the terms (including the 
source) in Eq. (0) by a factor 2p~^'S. After this rescaling, Eq. (^ takes the following familiar form 



fV-^") = 2p-4s{54[V''""^\9tV^"~^'] +r[matter]}. 



(11) 



In Ref. the wave operator was transformed to act on the field V'''^''=(p'"') '^'/'i'^^ rather than iIj^^'' (in order 
to achieve separability of the variables in the frequency domain) and takes the following form, in Boyer-Lindquist 
coordinates {t, r, 1!), ip) and Kinnersley tetrad 



r = 



(r^ + a'^f 



A 

"a(r - M) 



- c? sin^ d 
1 



a. 



\Mar 



dtp -4 



sin I? 

i cos-d 



1 



sin^i? 



a 
A 



sin^i? 



(4cot2i? + 2) 



1 cos i9 



(12) 



where M is the mass of the black hole, a its angular momentum per unit mass, E = + cos^ i?, and A = 
- 2Mr + a^. Note that if one wants to act on 

^/,(2) = (p(0))-4^(2) ^^^^^^ ^^^^ ^(2) ^^^^ 

one should 

consistently rescale all the terms (including the source) by a factor 2(p'^°))^''E (see Eq. (|ll|)). 

It is easy to show that a similar equation to ( pi] ) can be obtained for the Weyl scalar field V'O: upon exchange of the 
tetrad vectors I <-> n and m ^ m. In this paper we will explicitly work with -04 since it directly represents outgoing 
gravitational radiation. Since at every level of the hierarchy of perturbations we have the zeroth order wave operator 
acting on i/j^"-* we could always use the method of full separations of variables. In this paper, however, we will not 
proceed so because we want our equations to be suitable for evolution in the time domain from given Cauchy data. 



III. PRACTICAL ISSUES 



A. Gauge choice and computation of the source 



As we will show explicitly in the next Section, tp^ is neither invariant under first order coordinates transformations 
nor second order tetrad rotations. Thus, in order to integrate Eq. (|l^), one would have to evolve ■0 in a fixed gauge 
(and tetrad) and then compute physical quantities, like radiated energy and waveform, in an asymptotically flat gauge. 
This sort of approach was followed in Ref. to study second order perturbations of a Schwarzschild black hole in 
the Regge- Wheeler gauge which is a unique gauge that allows to invert expressions in terms of generic perturbations 
and thus recover the gauge invariance. There is not a generalization of the Regge- Wheeler gauge when studying 
perturbations of a Kerr hole, essentially because one cannot perform a simple multipole decomposition of the metric. 
Instead, Chrzanowski |32[ found two convenient gauges that allowed him to invert the metric perturbations in terms 
of potentials ^irg or '^org satisfying the same wave equations as the Weyl scalars p~'^ip4 or ipo respectively. 

In the ingoing radiation gauge (IRG) 

= = = = ^ = = = h^'L, (13) 

the homogeneous (for vacuum) metric components can be written, in the time domain, in terms of solutions to the 
wave equation for p^^?/;!^^ only, as follows 
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ihi'^) IRG = 2Re 



{~lt,lv{5 + a + 3/3 - r)((5 + 4/3 + 3r) - m^m^{D - p){D + 3p) 



+^(^rn,) [{D+^- p){5 + A(5 + iT) + {5-a + i(5-^-T){D + 3/9)] } {-^irg) 



(14) 



where Re stands for the real part of the whole object to ensure that the metric be real |3^^^J| and we made the e = 
choice. Note that in this gauge the metric potential has the property to be transverse (/i^^^^^ — 0) and traceless 
{h\P = 0) at the future horizon and past infinity. This is thus a suitable gauge to study gravitational radiation 
effects near the event horizon. 

The complementary (adjoint) gauge to the ingoing radiation gauge is the outgoing radiation gauge (ORG), 

/^ll^ - = /iiV = = h\ll = = h^^l = = /lil. (15) 

where the metric potential has now the property to be transverse (h'"^}n^^ = 0) and traceless {h'j}^ = 0) at the past 
horizon and future infinity. It is then an example of a suitable asymptotically flat gauge to directly compute radiated 
energy and momenta at infinity (see Sec. III.C). In this gauge, the homogeneous metric components can be written 
in terms of solutions to the wave equation for 4'q^\ as 



p ^ {-n^n^{5 - 3q; - /3 + 57r)((5 - 4q; + tt) - mf,m^{A + dp, - 3^ + j){A + p - 47) 



+n(tJri^) [{S - 3a + (3 + 5tt + t){A + p - Aj) + {A + 5p ~p - 3-/ - 'J){S - ia + tt)] } {^org) 



(16) 



Note that Eqs. ( p^ (or (|T5|)) are four conditions on the real part of the metric. Although (|l^)) (or (|l6|)) do 
not fix completely the gauge freedom, Chrzanowski metric choice given in Eq. ([l^)) (or Eq. (|l^)), being a specific 
choice between all the possible solutions satisfying those conditions, does uniquely fix all of the extra freedom. 

The potentials ibg and '^org fulfil the Teukolsky equation for p^'^ipi a-nd V'o respectively. To determine them we 
can invert expressions Eqs. (^3|) or (jl^) and its time derivatives at the initial Cauchy surface to relate the potential 
to our first order initial data. Alternatively, one can use the relations of these potentials to gauge invariant objects 
like ^po or p~'^'ip4. For instance, in the IRG we can take the relation -00 = DDDD"^ jug (See Eq. (5.28) of Ref. [|4|) 
or in the ORG the adjoint relation -04 = AAAA^'o/jg- Here we lower the order of the time derivatives of ^E* to first 
order ones by repeated use of the Teukolsky equation potentials fulfill (See, for instance, Eq. (5.20) of Ref. [Q). Since 
one can always make a mode decomposition of the ip dependence one ends up with a set of potential equations for 
^(r, 9) and 9t^'(r, 9) at the initial time. Boimdary conditions are chosen such that we obtain bounded solutions. The 
numerical integration of these equations is left for a forthcoming paper |4^ . These solutions give us the initial data 
to integrate the wave equations and then build up metric perturbations form Eqs. ( [l^ ) or (p^. The imposition of 
initial data to "04 and ipQ is discussed in the next subsection. 

Finally, in order to integrate Eq. ( ITTl ) we assumed the knowledge of the source term (^ since it depends only on 

first order perturbations. In practice, one solves the Teukolsky equation for (and / or ipQ^^) and builds up metric 
perturbations. It then remains the task of writing all first order Newman-Penrose quantities in terms of /i^^. This 
is not a trivial task, so we give all the equations relating the Newman-Penrose fields to the metric perturbations in 
appendix A. 



B. Imposition of initial data 



To start the evolution one has to be able to impose initial data to the second order invariant waveform. We first 
note that, from its definition, we can write 

^('^ - -c''^ + -h^^^h^'^ f ^("^ + - 2 ( h^^^ ~ -h^'-'^ \ ^('^ - 2/1^') ^^^^ im 

For the sake of definiteness we have used here Eq. ( |AlD choice of the first order tetrad, but it is clear that the above 

— . . . . . (2) 

expression can be written in a generic tetrad. Besides, since we are going to build up the invariant 0} , any choice 

of the tetrad (and the gauge) leads to the same, correct, result. 
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In Ref. 
expressior 



we have completely expressed V'i^'' (and its time derivative) in terms of hypersurface data only. The 



n'-m fi^m + 87V 



m n m 



- [0^1 - [0— '] 

w m w m 



(18) 



and its time derivative hold in general, to all order. Here N = (— g**) — N'^g*^, = n'' + N^n^ and 

TO = + N^m*. When we expand the above relation to a given perturbative order n, the proof given in Ref. p^ ] 
implies that ip4 and dtip4 will be independent on the lapse and shift of order n (but will depend, of course, on all 
lower perturbative orders of TV and N'^). 



To express our second order object 
hypersurface we will proceed as in Ref. 
perturbations. We then find 



/,(2) 



in terms of the three-metric and the extrinsic curvature of the initial 



26| taking now into account the additional terms, quadratic in the first order 



nmnm 



Rijkl + 2Ki[kKi]j 



TO h'^m + 8Af(o) 



(2) 



-4^(0) 
-87V(i) 



K 



■j[k,l] 



(1) 



7V m n m 



-87V(o)7V(i) 



. [0^1 - [0— '1 



(2) 



m n m 



. [0^1 - [0— '1 



(1) 



(19) 



(2) (2) 

Note that the first three terms have the same structure as in the first order case [for terms linear in h)/ and K}, .] 

L <'J ^ 

There is not dependence on the second order lapse and shift, but N(^i^ and the perturbative shift now explicitly appear. 
To re-express them in terms of hypersurface data, we can make use of Eq. (|l^ ) and Appendix A expressions that 
relate all first order quantities to directly expressible in terms of hypersurface data only as discussed before. 

And the same technique allow us to build up the additional quadratic terms occurring m. ip] '. Since the total ip\ ' 
was originally invariant, its final expression is not affected by the use of the a gauge choice (such as (14) or (|l^)) at 
an intermediate step. 

(2) (2) (2) 

For dtipj , the procedure is the same as before. We note that terms linear in hl^ and will have the same 



structure as in the first order case, so Eq. (3.2) of Ref. |26| apphes upon change of the subindcx (1) by (2). The 
additional terms, quadratic in the first order perturbations, can be directly written in terms of dt^ irg by taking the 
time derivative of Eq. (nlh and Appendix A expressions. 

In Appendix B we give an independent derivation relating to the four-geometry. We split 



= 



(20) 



where the first term on the right hand side is linear in the second order perturbations of the metric, i.e. hli^ and is 
formally the same as V'i^'' replacing h^i^} /i^V- The second term on the right hand side, i.e. "iAlq: accounts for the 
quadratic part in first order metric perturbations. 



C. Radiated energy and momenta 



The energy and momenta radiated at infinity to second perturbative order can be computed using the standard 



methods of linearized gravity (here h^i, stands for h 
null infinity). For outgoing waves [pO[ 



(1) 



defined in asymptotically flat coordinates at future 



^Note that the factor 8 appearing in front of the second bracket corrects an obvious misprint in Ref. 



This also applies to 



the Eq. (3.2) for dtip. 



(1) 
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lim V4 = ~7^{dth^^ 



(21) 



the total radiated energy per unit time (u ~ t 



dE 

—— = lim { - — 

du r- 



r) can thus be obtained from the Landau-Lifschitz pseudo tensor as 
du ip4(u,r,'d,ip) }, dn = sin -d dd d(p, (22) 



where we can consider ip^ — ijj^ 



(1) 



Note that Eq. (g2|) can be equivalently obtained by directly calculating the Bondi definition of the mass carried 
away by the gravitational radiation by imposing asymptotically flat conditions to the full Newman-Penrose quantities 
in a general vacuum spacetime. In this way, one can also compute the total linear momentum radiated at infinity per 
unit time along cartesian-like coordinates as [p5| 



dP^ 
du 




dn /„ 



du ipi{u, r, -d, if) 



(23) 



= (1, — sin 6* cos ~ smWsm(/3, — cosf 
and the angular momentum carried away by the waves |pq] can be obtained from 



dJ^ 
du 



lim 



47r 



Re 



dfi a 



du V-'4(w, r, d, If) 




du ip^iu, r, d, (f) 



(24) 



One can directly compute the second order correction to the energy and momentum radiated at J'^ using V'4^\ 
provided one is working (to first order) in an asymptotically fiat gauge (for instance, the outgoing radiation gauge). 
Eqs. (^2|)-(|2j), written in terms of the full, nonlinear ^/'4, are covariant expressions, holding in any asymptotically flat 

spacetime. To first perturbative order, ip'^^^ is directly gauge and tetrad invariant, so one can forget that the above 
equations had been obtained in an asymptotically flat gauge and think of them as gauge (and tetrad) invariant. We 

(2) 

would like to have the same nice property to second perturbative order, but ^\ is not invariant. One should then 
build up a gauge (and tetrad) invariant waveform -0^^'' that, in an asymptotically flat gauge coincides with -0^^^ 
This will ensure us the direct use of Eqs. (]2^)~(p^) in terms of our invariant object, i.e. 0^^^ given in Eq. (^). 



IV. CONSTRUCTION OF THE SECOND ORDER COORDINATE AND TETRAD INVARIANT 

WAVEFORM 

The general covariance (i.e. diffeomorphism invariance) of Einstein's theory of gravity guarantees the complete 
freedom in the choice of the spacetime coordinates (gauge) to describe physical phenomena. In the relativistic 
theory of perturbations one always introduces two spacetimes, the physical (perturbed) spacetime and an idealized 
(unperturbed) background. In this way the perturbations can be viewed as fields propagating on the background. 
Consequently, to compare any physical quantity in the perturbed spacetime with the same quantity in the unperturbed 
spacetime it is necessary to introduce a diffeomorphism about the pairwise identification points between the two 
manifolds. The arbitrariness in the choice of this point identification map introduces an additional freedom to the 
usual gauge freedom of general relativity and is at the origin of the gauge problem in perturbation theory ]37| . A 
convenient way to deal with this gauge problem is to construct quantities which are invariant under a change of the 
identification map of the perturbed spacetime while the background coordinates are held fixed. 

Invariance in the Newman-Penrose formalism has a more restrictive meaning than in the standard (metric) pertur- 
bation theory, since the introduction of a tetrad frame at every point of the spacetime now requires that any physical 
perturbation must be invariant not only under infinitesimal gauge transformations (GI), but also under infinitesimal 
rotations of the local tetrad frame (TI). In this Section we briefly review the basic concepts of (higher order) tetrad 
invariance and coordinate (gauge) invariance in the framework of the Newman-Penrose formalism. We start with our 

second order object ^\ , which is neither invariant under first order changes of the coordinates nor under second 

(2) 

order tetrad rotations. We then show how to build up a tetrad invariant object by adding to ^4 a conveniently chosen 
term, quadratic in the first order perturbations. In this way the new object will be invariant under the (6-parameter) 
tetrad rotations. The procedure for the construction of the totally invariant object, i.e. also under coordinate choices 
(4-parameters) is analogous, but algebraically more involved. The final result is a general prescription for constructing 
totally invariant (I) quantities directly related to the (outgoing) gravitational radiation. 
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A. Tetrad invariance 



The 6-parameter group of homogeneous Lorentz transformations, which preserves the tetrad orthogonahty relations 
If^ = —ra^m^ — 1 (and aU other scalar products zero), can be decomposed into three Abelian subgroups: 



• Null rotation of type (I) which leaves the unchanged: 

Z^ ^ Z/i 



+ a,m^ + a.m^ + a.al^, (25) 
fa,, — > m,, + a/,,; 



Null rotation of type (II), which leaves the unchanged: 



n^, (26) 



Boost and rotation of type (III): 



1-1, 



^ A-'n^, (27) 
TO^ exp(z6')m^; 

where (a, b) are two complex functions and (A, 9) two real functions on the four dimensional manifold, hence the six 
arbitrary parameters. When these functions are taken to be infinitesimally small the above transformations can be 
expanded up to an arbitrary order and then applied to any Newman-Penrose quantity. 
Under a combined tetrad rotation of classes I, II, and III 

^ 4'^ + 2[(A - 1) - iO]iji^^ + 4a?A^'^ + 6aVr • (28) 

The idea here is to supplement t/j^ with additional terms that make the whole object tetrad invariant. Since we 
have to add those "correcting" terms on both sides of the field equation (|l^), we will write them as powers of first 
order perturbations so they can be added to the source term The first step towards constructing this quantity is 
to note that {Fy{m^)'^^4/[{F '^*'^)^(m'' (''))^] is invariant under rotations of class III. The second order piece of this 
combination of fields is 

(2) 

Note that the second addend exactly compensates for the variation of class III of tfAJ (proportional to the parameters 
A—1 and 6 in Eq. (p8|)). In addition one can easily check that the second term in (|2^) is also invariant under rotations 
of class I and II with the Kinnersley choice of the zeroth order tetrad 

(to^)(0) = — (iasin??,0,l,i/sin??) , 

V2(r + mcos?9) 

since F m"^ and ip^l'^ are all invariant^ under rotations of class I and II. Still, the first term in ( p9| ) varies with 
respect to rotations of class I and II. To correct that we note that under combined rotations I, II, and III 



^It is clear that we can write the tetrad invariant object in terms of a generic zeroth order tetrad by replacing in Eq. (32) 

ir (1) ^ [r (1) _ m'"V'i^V(3?2''') - m''i/>$^V(34°^) and - Z''V'i^V(3V'2°')- We take the background tetrad @ 

for the sake of simplicity. 
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^(i)^^W+3aV^f . (31) 

This allows us to solve for a and replace it into the new expression (its form suggested by the a dependence in the 
transformation (|2^)) that supplement (p9|). [Note that this replacement is successful because ip^ vanishes to zeroth 
order.] Thus, the object, 



is second order tetrad invariant. 

While the above combination is tetrad invariant, one can see from the general behavior of the Weyl scalars and spin 

coefficients in an asymptotically flat gauge (see, for instance Sec. VII of Ref. [0), that the quadratic term we added 

(2) 

does not vanish relative to ip\ for large r, i.e. goes like 0{l/r) as well. In order to have the desired property that 

in the radiation zone the invariant object approaches ipl^ {AF stands for an asymptotically flat gauge), we will 
subtract to ( |3^ ) another quadratic part that both, cancels its added asymptotic behavior and is tetrad (and gauge) 
invariant in order to preserve the gained invariance of (|3^). Symbolically, if we call Q the quadratic part we added to 
-04 in Eq. (p2|), we search for a 

0(%j=^f) +Q-Qf^ (33) 



A practical way to build up Qf^ is to use relations (|l6|), i.e. the perturbed metric in the outgoing radiation gauge, 
which is an asymptotically flat gauge at infinity. In this gauge, we evaluate the quadratic part Q in (^ and once, 



re-expressed all in terms of VPoflG via Eqs. (16), we can forget that we used the outgoing radiation gauge and see 



qORG tetrad and gauge invariant object, since ^org is totally invariant. In the outgoing radiation gauge ij}^^^ 
(and i/"! Ti) reduces to i/'l l ^ ^e directly deduced from the expressions given in Appendix B. 



B. Gauge invariance 



The meaning of gauge invariance under infinitesimal coordinate changes, to an arbitrary order in the perturbations, 
was explicitly elucidated in Ref. |Q following the approach of Ref. Q. Locally, these gauge transformations are 
the 4-parameter group of the inhomogeneous Lorentz transformations. Up to second order in the perturbations an 
infinitesimal change of coordinates 

where ^'^^^ and ^^2) ^^e two independent arbitrary vector fields and e a small (perturbative) parameter, produces the 
following effect on the first and second order perturbations of any quantity $ (scalar, vector or tensor field) that we 
assume can be expanded as ^^^^ + ^'^^^ + + ... 



$(1) ^ $(1) +£^^^^$(0)^ (35) 




where, for the sake of completeness we recall here explicitly the basic coordinate expressions of the Lie derivative 
along a vector field 

£^<I> = , if $ is a scalar; 

A*'' = ^r^^^ ~ C*^ ' if i« a sector- (37) 

C^^al3 = $a/3,M^'' + tJ.f3 + ^^^a^ ' *a/3 is & tenSOr. 

Note that from transformation ( ^5| ) it follows that all Newman-Penrose quantities that vanishes on the background 
(or more precisely satisfy /Cj^jj^'"^ = 0) , like f/'o^'') i^']l'\'4>'^\k^^\a^^\X'^'^\v^^\ are first order gauge invariant 
(GI). Transformation (^), however, states that none of these Newman-Penrose quantities, to the second order in the 
perturbations, are gauge invariant, since /Cj^^j*'-^' 7^ 0. Thus, none of the interesting Newman-Penrose quantities 
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that are tetrad invariant (TI) and gauge invariant to the first order are also invariant to the second order in the 
perturbations. In particular, second order gauge invariance requires that the quantity vanishes to zeroth and to first 
perturbative order. 

Explicitly, for the scalar field 1^4 we have 



',(0) „;,(2) 



Hence we see that the vanishing of tp^ ensures that t/j^ will be gauge invariant under "pure" second order changes 

of coordinates, but since ip^^"^ will in general depend on all four coordinates, tp^^ will not be gauge invariant under 
first order changes of the coordinates. 

In order to apply similar techniques to those we used to construct a tetrad invariant object now in the coordinates 

(2) 

context, i. e. by "correcting" ip\ with products of first order quantities, we will make use of the following Lemma 

Lemma: The product of the first order pieces T'-^^P^^^ of two tensors (that can be expanded into perturbations) 
transforms under a first plus second order gauge change, given by Eq. (Qj, as the product of the first order transformed 
quantities individually. 



(r(i)p(i)) ^ (r(i) + £5,^,r(°)) (pw + 



Proof: Let T and P be two general tensor fields. Apply the first plus second order transformation ( |36| ) to the 
product and consider second order pieces, then 



(TP) ^ (TP)(2) +£^^^^(rp)(i) + -{cl^^ +/:e^,,)(rp)(«), 



or more explicitly 

('7-(2)p(0) T(r)p{l) ^ p(0) p(2)) _^ ^p(2)p(0) _^ p(l)p(l) ^ p(0) p(2)-j ^ ^rj,{l) p(0) ^ j.(Q) p(l)-| 

)(T(")P(")). 

We now apply the same transformation (^) to the products P'^'^^T and T^^^P to obtain 
p(5^(2) ^ p(0)T(2) +p(")£^^^^(T(i)) + lp(o\cl^^ +£5,,,)(T)(°) 

similarly 

T{^{2) ^ r(")p(2) + r(")/:^,,,(p(i)) + ^t^"\cI^^ +/:^(,,)(p)(°). 

Upon subtraction of the last two expressions from the first one, we obtain 

TWi) ^ T(i)pW + T(i)£^(^j (P(°)) + P('^/:C(„ (T^")) + £5^^, {T^°^)C^^rj (P'^^^)- (39) 

This proves our Lemma. An obvious corollary is the case when both fields are gauge invariant, i.e. C^^-^^{T'^^'>) — 
and Z^^j^j (P^")) = this generates a second order quantity that is first and second order gauge (coordinate) invariant. 

(2) 

To construct a second order gauge invariant waveform ip^ qj we can then use the same techniques as in the previous 

subsection. It is convenient now to start from our tetrad invariant object, as defined in Eq. (pq). Under a first order 

(2) — 
coordinates change ip^ rpj transform as 

V4 TI -*VaTI+ + ^^4 I ir (0) + (0) I ' 



where we made use of the properties expressed in Eqs (pq) and (|39|) 
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As in Section III. A, the idea here is to add to ^'''Pti terms quadratic in the first order perturbations in order to 
make the whole object coordinate invariant^ while preserving its tetrad invariance. The procedure can be summarized 
as follows, 

Prescription: The first step is to invert the coordinate transformations of first order quantities for the gauge 
vectors ^^^j^^. We shall denote this first order combination by the boldface vector: ^^ly i-e. ^^^^ = ^^^-j — Cj^j^-). Making 
the replacement ^^j^^ ^^(i) ( F^ above generates a totally invariant object. Still from all the possible 

invariant objects we want those whose quadratic term do not contribute to the radiation in an asymptotically flat 
gauge (AF). As we discussed at the end of Sec. IV. A, this ensures us a simple interpretation of the invariant -0/ 
regarding radiated energy and waveforms. Since Eq. ( pO| ) is linear in Sj^jy subtracting the quadratic term in an 

asymptotically flat gauge will be equivalent to make the following replacement ^^^^^ ^(i)*^ ^ ^(i) ' discussed 
before, a practical way to evaluate ^^^^'^ and keep the tetrad and coordinate invariance is to use the outgoing radiation 
gauge (Eq. (p^) and consider the final expression in terms of "^org a-s a totally invariant expression regardless its 
derivation with a choice of the first order gauge and tetrad. 

We recall here that '04 and of course also terms quadratic in the first order perturbations are already invariant 
under pure second order coordinate transformations, labeled by ^^j) • 
symbolically expressed as 



Finally, our invariant waveform can then be 



ORG 
^(1) 



r,r (0) 

I', LI 



( ORG A Iti (0) ( fr ORG 



(0) 



+- 



ORG 



- 

^(1 



ir (0) 



ORG ^-d 



(0) 



(41) 



C. Construction of the second order invariant waveform 



The above prescription is conceptually very simple. However, in practice, to find and €''(1) brings some 

technical complications. The first remark is that the procedure is not unique. We have a big choice of first order 
objects (all Newman-Penrose quantities, metric, extrinsic curvature, etc) to build up ^'(^ly In fact, one can easily see 
that the ambiguity to generate an invariant waveform has to be present since one can always add products of first 
order invariant objects to generate a new second order invariant object. The requirement that the quadratic correction 
must not influence the asymptotic behavior greatly reduces this ambiguity. In fact, physical quantities such as the 
radiated energy and observed waveform, defined in an asymptotically flat region, are uniquely defined by this method, 
since the differences introduced by different asymptotically flat coordinates vanish with a higher power of r. We thus, 
give an explicit object in order to be able to make comparisons with, for instance, full numerical results that directly 
compute the covariant object ip4. Below we give a simple choice of ^'^ly in order to construct that is valid for 
perturbations of Kerr black holes, i.e. a 7^ 0. In Appendix C we give another choice for the case of a Schwarzschild 
background. 

In the rest of this subsection, to simplify the notation, we drop the subscript (1) from the first order gauge vectors 
since we will never refer to the second order gauge vectors. The and ^ components can be easily found from 
the variations of the tetrad invariant Weyl scalar 



(1) 



(1) 



(1) 



(42) 



and of its complex conjugate 02 



(1) 



c = 

^ 6M 



02 



0(^) 



(43) 



*A similar procedure was adopted to generate second order gauge invariants in the Moncrief s formulation of Schwarzscliild 
black hole perturbations Un. 
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6M(msini^) 



V2 



4'' 



(44) 



The same techniques cannot be straightforwardly appUed to find the other two components ^* and ^'^ . The origin 
of the problem can be traced back to the fact that the Kerr metric has two two killing vectors along dt and d^, 
and thus one can never find local, first order quantities that vary with ^* or but only with the derivatives of 
them. Explicitly, using the variations of the metric and extrinsic curvature components (which are tetrad invariant 
quantities), we find (here background fields are unlabeled) 



9w 9rr^,t) 9tip{hrip 9rr^,ip) 

2 ' 

9t:p ^ 9tt9vv 
{htil + 9i)^i^t + 9tvCi,) 
9tt 



(45) 



and 



gtt{hi]) + grrCJ - gtvihil^ + grrCt) 



K. 



(1) 



9tip ~ 9ttg(p(p 



(46) 



Thus, to find ^* and one has to integrate their four derivatives over the spacetime. This can be performed like 
the integration of a potential in four dimensions. For that one has to verify the integrability conditions. In practice, 
since we are going to compute differences of these vectors, with respect to the asymptotically flat ones, the existence 
of the ^* and components are assumed a priori and they are related to the existence of the outgoing radiation 
gauge proved in Ref. |Q. As a consequence of these integrals on first order fields, the resulting waveform will be 
nonlocal, but this carries no further consequences since in solving the second order perturbations we assumed first 
order ones to be completely known. Notably, the evolution equation for the second order perturbations is local. In 
fact, only derivatives of ^* and enter in building up the source 



T[4P] - Si 



^f)-(p(«))-vn, 



where the source term (as can be derived from Eq. (^) is now 
Si = 2(p("))-4s {Si[4l^] + ^[matter]} 



(47) 



(48) 



-r 



2i) 



(1) 



/r (1) _ /'■ (1) j„i? (1) _ (1 

ORG li^nur 



ORG 



ir (0) 



7,^ (0) 



^(1) ORG 



(p(o))Vr 



r 



^,fj. ^ (1) 



(1) 



''A* S (1) 



C (1)^" 



(1) 



ir (0) 



m 



■d (0) 



(49) 



here C^^^f^f - If^^. 



(!)• 
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While the evolution is local, we need to compute the waveform at least on the initial hypersurface and then at 
the observer location (to compute, for instance, the radiated energy). At t — 0, after mode decomposition in the 
coordinate, we have 



/^i i-i ORG\ 

^'^■V ^-V >™ _i / ^^(ti ti ORG\ I / ^^qfci ORG\ / ^^^,qfti ti ORG\ 

^-^ zm 



where i = {t, tp) and the same equation holds for the observer at a fixed rots, exchanging the roles of r and t. 

Note the presence of the integration constants c* and . They represent first order changes in the origin of time and 
azimuthal angle. This problem was already found in Ref. |^ and there it was given a method to fix "a posteriori" the 
value of the constants. In Section IV. B we generalize the procedure given in Ref. ||^ and cxplicity write the integrals 
that are necessary to fix the constants c* and c"^. 

(2) 

In order to compute the totally invariant second order waveform we must fix the constants ct and c<^. We can 
generalize the gauge fixing prescription given in Ref. and define "a posteriori" ' the value of the constants 

C^dt{4'^rC^dt{d,4'^r-j-^dti.i'^dA^^ 



C^dt{^.^y)^ s°^^dt{d,^^^Y ^ j-'^dt^^^ 

We can then construct the "c-invariant" waveform 

i^fc = V-f - {c'dt + c^d^)^^^^ , (53) 

This procedure amounts to gauge fixing the zero of time and of the azimuthal angle in such a way that the integrals 
in the numerators of Eqs. (|5l| ) and ( p2[ ) vanish. In order to be able to compare the perturbative results with the full 
numerically ones it is crucial that one is able to perform the same origin of coordinates fixing. Note that we can also 
fix these constants at the initial hypersurface t = Q. The same expressions (|5T|)-(|5^ apply changing the integrations 
in time by integrations in r. 



V. SUMMARY AND DISCUSSION 



In this paper we presented a gauge and tetrad invariant framework for studying the evolution of general second 
order perturbations about a rotating black hole. To do so, we first uncoupled second (and higher) order perturbations 
of Kerr black holes for the Weyl scalar '04, that directly represents the outgoing gravitational radiation, and found 
that the perturbed outgoing radiation field V'i"'' fulfils a single Teukolsky-like equation (see Eq. ([Ill)) with the same 
wave operator as for the first order perturbations [ pO[ , acting on the left hand side and an additional source term 

— ... (2) 

written as products of lower order perturbations on the right hand side of the equation. We note, however, that ipX 
is neither tetrad nor first order coordinate (gauge) invariant. It is only invariant under purely second order changes 
of coordinates, simply because 1/^4 vanishes on the background (Kerr metric). Invariant objects to describe second 
perturbations lead us to reliable physical answers without having to face gauge difficulties. Hence, we explicitly 
show that it is always possible to correct '04 in order to build up a complete second order invariant waveform '0^ 
(i. e. invariant under both tetrad rotations and infinitesimal coordinates transformations) that gives a measure of 
the outgoing gravitational radiation. This is done in Sec. IV where we give a general prescription to produce the 
result expressed in Eq. (ETl). We also show that the same equation as (O), with a "corrected" source term, is now 
satisfied by (see Eq. (E^). A number of interesting conceptual and technical issues raised from this computation, 
like the appearance of nonlocalities in the definition of the gauge invariant waveform when we want to relate it to 
known first order objects and its non uniqueness. Seen in retrospective, our method of generating a gauge invariant 
object is like a machine that transforms any (first order) gauge into an asymptotically fiat one, in particular, into 

the outgoing radiation gauge. In fact, in this gauge we have (04^^)'^''^'^ = (0'4i)'^'^'^ = {ip^P)'^^'^ = ■0/^''. This fits 
into Bardeen's interpretation of a gauge independent quantity and suggest to work in the outgoing radiation 
gauge as a particularly simple way of dealing with the numerical integration of the second order equations p2| . In 
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the language of Eq. (BO) we see that the process of building up V'/ is like subtracting the first order piece to 'tp^. Our 
gauge invariant object, V'/, is not the second order term of a series expansion of "04, but it can be related to psi\ in 
an asymptotically flat gauge. 

The spirit of this work has been to show that there exists a gauge invariant way to deal with second order per- 
turbations in the more general case of a rotating black hole and to provide theoretical support to the numerical 
integration of the second order perturbation problem. In order to implement such integration of Eq. ( ^7| ) we proceed 
as follows: We assume that on an initial hypersurface we know the first and second order perturbed metric and 
extrinsic curvature. We then solve the first order problem, i.e. solve the standard Teukolsky equation for ■tp^'' (and 

for V'o^^). Next we build up the perturbed metric coefficients in, for instance, the outgoing radiation gauge ( p^ . The 
perturbed spin coefficients are now given by expression (A4) and the perturbed covariant basis by (A2). Those are 
all the necessary elements to build up the effective source term appearing on the right hand side of our evolution 
equation, as explained in Section III. A. It is worth to note here that from the analysis of the asymptotic behavior of 
the different Newman-Penrose quantities |l3| involved in the source, one can see that at infinity the envelope of (the 
oscillating) S4 is at least of 0{r~'^), which guarantees the convergence of the integration of Eq. (p7|). 

The other piece of information that we need in order to integrate Eq. ( p7| ) is i/'j^-' on the initial hypersurface. This 
is explained in Section III.B. We also need to use in this case Eq. ( po|) and the expressions given in Appendix B. For 
the computation of the radiated energy and momentum one uses Eqs. ( p2| ) and (|2^). The advantage of this procedure 
is that we can now use the same (2+l)-dimcnsional code for evolving the first order perturbations by adding a 
source term. In fact, the background (Kerr) metric allows a decomposition into axial modes, i.e. the variable ip. A 
mode decomposition of all quantities involved in the second order evolution equation can be trivially performed (note 
that in the source, involving quadratic terms in the first order perturbations one has to include a double sum over 
modes, let us say, m and m'). In the time domain no further mode decomposition (i.e. in multipole) of the source 
term is practical. 

An important application of the formalism presented in this paper is to reproduce the results obtained in Refs. 
(for the nonrotating case and the multipole £ = 2). The complexity of the calculations in the standard Zerilli 
formalism that would follow from considering the sum over all multipoles can be notably simplified in the Newman- 
Penrose formalism. We can thus also study the £ = 4 multipole of the radiation and not only test the efficiency of 
our formalism, but also make a more detailed comparison with full numerical results. The next step is to extend the 
numerical computation to the more interesting case of rotating black holes. The numerical integration of Eq. ( ^7| ) will 
be relevant not only for establishing the range of validity of the collision parameters in the close limit approximation, 
but (hopefully) to produce a more precise computation of the gravitational radiation. Direct comparison with the 
existing codes for numerical integration of the full nonlinear Einstein equations is possible . 

Following the steps described in this paper, upon exchange of the null directions / <-> n and m <-> m, it is straight- 
forward to write the corresponding equations for ipQ in case one wants to have a description in terms of ingoing 
waves. This would allow to study the influence of gravitational radiation on the horizon of a rotating black hole, 
critical collapse and also phenomena in their interior, like the mass-inflation. We studied in detail only gravitational 
perturbations, but it seems straightforward to generalize our method to scalar and vector perturbations. We also note 
that, although we have focused our attention on the problem of colliding black holes, the second order perturbative 
formalism developed in this paper can be easily generalized to any Petrov type D (or even type II) background metric 
and thus can be applied to study other interesting astrophysical scenarios as nonrotating neutron stars and cosmology. 
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APPENDIX A: FIRST ORDER NEWMAN-PENROSE QUANTITIES 

Throughout this Appendix, to simplify the notation, we omit the superscript (0) on the background quantities, while 
all first order quantities are denoted with the superscript (1) with the exception of the first order metric perturbation 
that we simply denote as /i^j/. 

Let us first note that the perturbed null tetrad can be represented by [Q 
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;f(i) ^ --hun'', 
2 



,M(i) _ 



(Al) 



h 



Note that in order to have this cxphcit form a choice of the first order null directions was made. To relate this to the 
metric perturbation recall that 5^1/ — 21(^^71^^ — 2m(^(TOi,) which implies that /i^i^ = 2Z|^-'n^') + 2l(^^rS^^ — 2r7^|^■'m^,) — 



Making use of the relations (Al) we can immediately derive the first order Newman-Penrose directional derivatives 



^(1) 


= = - 




Ad) 






^(1) 







(A2) 



In order to compute the spin coefficients to the required order we follow Ref. 
relations ||3l[], Ch. 1.8 (these are exact expressions) 



making use of the commutation 



AL* - DA = (7 + 7) L> + (e + e) A - (t + tt) (5 - (t + 7f) ^, 
5D - D5 = {a + (3 -Ti) D + kI^ - {-p + e - t) 5 - aS, 
(5A - A(5 = -VD + (t - a - /3) A + (/I + 7 - 7) (5 ~ A5, 
M - 55 ^ {T^ - D + {-p - p) + {a -13) 5 + {P ~a) J, 



(A3) 



expanding both sides to first perturbative order and using Eq. (A2) we can equate the coefficients of each operator 
to get a system of linear equations (16 of which only 12 are independent) that can be solved for the spin coefficients 
giving 



ad 

Ad 

2p^- 

26d 
2^d 

2rd 
2a(i 



(D-p-2 e) h 



1 



Im 



(^- 2a- 2/3 + 7f + r) hi, 



(tt + t) + - (D + p - p + 2e - 2e) 

- (A + p + 27) + i + 2a + 2/3 - TT - r) /i„„, 



(t + Tt) hr, 



- (A + 71 - /i + 27 - 27) hr, 



phnn - ((5 + 2/3 + r) ha— I + (5 + 2/3 - 27r - r) /(,„m - - (2 A + 71-/1 + 7-7) h^rcm 
JIhu + {p - p) hni + {D + p - p) hmm - {5 - 2a ~ If) hm + (J + 2r - 2a + tt) him, 
(D + p -p) hni + - (J - 2a - tt) hijn - ^ ((5 - 2a + Stt + 4t) hmi 

+ \ (P - P) hmrn - ^ (A + 27) hu, 



-{D - p-2e) hnrn - ((5 + T + Tt) /l„; - (A + /i - 27) hm ~ ThmWi - Thj^ 
(D -p+ 2e) hnm + (5 - TT - r) /l„( + (A + p - 27) him - Trhmrn ~ T^hmm 

^{D~2p- p-2e)hnm-^{A + 'ij-2p + -p- 27) hm 
-- + TT + r) + i + 2a - TT - r) /imm - - ((5 - 2a + 7f + r) /imm 
2^d) = i (D - p - 4e + 2p + 2e) - ^ (A + /i + 271 + 27) /i,™ 

((5 + 7f + r) ft,„( - i ((5 - 2/3 + TT + t) /i^ttj + i (^ + 2/3 - tt - t) 



(A4) 



16 



27(1) ^ - + hni + ^{D + p -p + 2e) hnn -^{S + 2p + 2W+ 3r) h^j^ 

+ i (5 + 2/3 - 27r - r) /i„„i + i (Sp - 2/i + 7 - 7) /i^^ji 

Note that these expressions are completely independent of the choice of the gauge, although a tetrad choice to first 
order had to be made in Eq. ( |Al| ). 
Finally, the exact Weyl scalars are 



■00 = (-D — 3e + e — p — p)a — ((5 — a — 3/3 + 7r — t)k, 



(A5) 



V'l = (-D + e - p)/? ~ {6 — a + n)e - (a + 7r)cr + (7 + p)k, 

■02 = [(J - 2a + /3 - TT - t)/3 - (J - a + 7f + r)a 

+ {D + e + e + p - p)-f — (A — 7 — 7 + 71- /i)e 

+ (J - a + ^ - r - 7r)T - (A - 7 - 7 + 71 - + 2{i^k - Act)] /3, 

ips = (S + ]3 -t)-/ - {A - ^ + -p) a + {e + p)iy ~ (P + t)A, 



and 



^4 =: (,5 + 3a + /3 + TT - T)f^ - (A - 7 + 37 + M + 7i)A, 



(A6) 



Note that these expressions can be trivially expanded to first perturbative order and hold when matter sources are 
included. 

APPENDIX B: SECOND ORDER NEWMAN-PENROSE QUANTITIES 

Taking the same first order choice of the tetrad to second order (we can do this because the final aim is to plug this 
into an invariant object) one obtains 



;m(2) 



,(1)^(1) 



If" 



i (2) 



2 nn ^ ^"'nm ''nm 



'Hi "'nn ^ 2 \ lm"'nm ^ "'lm"'nm 



(Bl) 



"nin ' ■"nm mm ' mm nm 



mm ' 2 mm mm ' 2 mm mm 



mm mm mm 



We now expand up to second order the third commutator of Eq. ([A^) to obtain z/^^^ ^ + ^ ^-nd A'^^^ ^'^l + 



i^i^) = -(A + 71 + 27)/ii'l + i(<5 + TT - r)/.!^^, 

= -(A + 71 + 27)/i« + 2(5 + ^ - t)/.« /^i'^ - A(i)/.W 



/3)^'^ + (7 + 7)^ 



Im 



(^_7 + 7)(i)_(7 + 7)/,W 



,(1) 



Af = -i(A + 71 - /i - 27 + 27)/iEr - + 



A 



(2) 



1 



(A + 71 - M - 27 + 2j)h!^hl^l + (r + tt) 



"'mn'nn "rnm"'nm 2 ml"-nn 



(B2) 
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Finally, using Eq. (^) we find t/.^ ^ = vIl + ^^^^ 



4Q' 

^('l' = + 3a + ^ + tt - r)4'^ - (A + 71 + /i - 7 + 87) , 

V;g = + 3a + + TT - T)i.J^ - (A + 71 + /i - 7 + 37)Ag^ (B3) 
+ (^+3a + ;5 + 7r-r)(iV(i) - (A + 7I + ^ - 7 + 37)(^U(^^ 



APPENDIX C: GAUGE INVARIANTS IN THE SCHWARZSCHILD LIMIT 

In the case when a = 0, we can find the following set of first order gauge vectors assuming that the is given by 
Eq. (H) 

cos-&{2h^^ + 2 sin^ "^r^l^ — sin^ + sm'd{2imh^^ — h^^^§ — 2sin^ 'i^'rCm^s) 



i = -- 



2r2 sin'!9(m2 + 1) 



" 2(r-3M) '^^^*'^''' ~ '^^^^ '^^ ^"^^ ^ ~ 2M)|*„ 



wTilJJ, = - ^ ^ . 2 a, 9 , -,^ {("1^ + 1 - 2cos2 +cos2i9sin2i?/ii)tf 



1 

2r2 sin^z9(m2 + 1) 

+ cosd sin ^{hipip^^ — 2imh^^ + sin^ 'J'^^m.ij) + 2 sin^ i9(m^ + sin^ i?)€m}7 



^m,^ = ~ sin 2^(r - 2M)(r - 3M) {"'"^'^'^ty-^ + m'^'^m,r " WMr^^^, - imr^ sin^ 



2imMrht^^r + Am^M'^^l^^ + 2imMR^ sin^ + 2miM/it^ - 6Mr sin^ ^Cm,t + 2r-^ sin^ 'd^';^^t 

sin 2-d{2,Mht§ - rhti)) - SMh^^^t + rh^^^tr} (CI) 



„ = - - — — — ^— - 2 i sin ?9 cos^ i^/i^j^ + 2 m cos ??/i^<p + 2 m cos sin^ t^?'! m 

(2 (m^ + 1) sm v) \ 

— 2 ih{)ip siv? i9 — m sin^ d cos dh^ — m sin dh^^^^ — 2m sin"^ ^^''^ m.i?^ 1 
2Afr2^:„ + htt 



2r(r - 2M) 



(r - 2M) [r/.,^,. + sir? ^^rt + Mr - ^C^J - ^Mh^^ 

2r2 sin^ i?(r - 3M) ^ 

^'^ and as given by Eqs. (^) and are well defined in the a ^ limit. 
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